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Abstract 

Eigenvalue problems are fundamental to mathematics and science. We present a simple 
algorithm for determining eigenvalues and eigenfunctions of the Laplace-Beltrami operator 
on rather general curved surfaces. Our algorithm, which is based on the Closest Point 
Method, relies on an embedding of the surface in a higher-dimensional space, where standard 
Cartesian finite difference and interpolation schemes can be easily applied. We show that 
there is a one-to-one correspondence between a problem defined in the embedding space and 
the original surface problem. For open surfaces, we present a simple way to impose Dirichlet 
and Neumann boundary conditions while maintaining second-order accuracy. Convergence 
studies and a series of examples demonstrate the effectiveness and generality of our approach. 
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1. Introduction 

The study of eigenvalues and eigenfunctions of the Laplacian operator has long been a 
subject of interest in mathematics, physics, engineering, computer science and other dis- 
ciplines. Of considerable importance is the case where the underlying domain is a curved 
surface, S ) in which case the problem becomes one of finding eigenvalues and eigenfunctions 
of the corresponding Laplace-Beltrami operator 

- Ys • V s u = \u, (1) 
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or, more generally, the elliptic operator 

— Vs • (a(x)Vsu) = Xu. 

The Laplace-Beltrami eigenvalue problem has played a prominent role in recent years in 
data analysis. For example, in [1], eigenvalues of the Laplace-Beltrami operator were used 
to extract "fingerprints" which characterize surfaces and solid objects. In [2, 3], Laplace- 
Beltrami eigenvalues and eigenfunctions were used for dimensionality reduction and data 
representation. Other application areas include smoothing of surfaces [4] and the segmen- 
tation and registration of shape [5] . 

Analytical solutions to the Laplace-Beltrami eigenvalue problem are rarely available, so 
it is crucial to be able to numerically approximate them in an accurate and efficient manner. 
Partial differential equations on surfaces, including eigenvalue problems, have traditionally 
been approximated using either (a) discretizations based on a parameterization of the surface 
[6] , (b) finite element discretizations on a triangulation of the surface [7] , or (c) embedding 
techniques which solve some embedding PDE in a small region near the surface [8] (see also 
the related works [9, 10, 11, 12, 13, 14, 15]). 

Parameterization methods (a) are often effective for simple surfaces [6], but for more 
complicated geometries have the deficiency of introducing distortions and singularities into 
the method through the parameterization [16]. Approaches based on the finite element 
method can be deceptively difficult to implement; as described in [7], "even though this 
method seems to be very simple, it is quite tricky to implement" . Embedding methods (c) 
have gained a considerable following because they permit PDEs on surfaces to be solving 
using standard finite differences. 

This paper proposes a simple and effective embedding method for the Laplace-Beltrami 
eigenvalue problem based on the Closest Point Method. The Closest Point Method is a recent 
embedding method that has been used to compute the numerical solution to a variety of 
partial differential equations [17, 18, 19, 20], including in-surface heat flow, reaction-diffusion 
equations, and higher-order motions involving biharmonic and "surface diffusion" terms. 
Unlike traditional embedding methods, which are built around level set representatives of 
the surface, the Closest Point Method is built around a closest point representation of the 
surface. This allows for general smooth surfaces with boundaries and does not require 
the surface to have an inside/outside [17]. In addition, the method does not introduce 
artificial boundary conditions at the edge of the computation band. Such artificial boundary 
conditions typically lead to low-order accuracy [12]. 

Here we apply the Closest Point Method to the problem of determining the eigenvalues 
and eigenmodes of the Laplace-Beltrami operator on a surface. We begin by demonstrating 
that, for closed surfaces, there is a one-to-one correspondence between the eigenvalues of 
the embedding problem and the original surface problem. Later, we consider open surfaces 
and present simple techniques for achieving high-order accurate approximations to Dirichlet 
and homogeneous Neumann boundary conditions. Our proposed method retains the usual 
advantages of the Closest Point Method, namely generality with respect to the surface, 
high-order accuracy and simplicity of implementation. 
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The paper unfolds as follows. Section 2 provides key background on the Closest Point 
Method. Section 3 proposes an embedding problem used to solve the Laplace-Beltrami 
eigenvalue problem and explains why a similar embedding problem leads to spurious eigen- 
values. Section 4 provides discretization details. In Section 5, a second-order discretization 
of boundary conditions is described for open surfaces. Section 6 validates the method with 
a number of convergence studies and examples on complex shapes. Finally, Section 7 gives 
a summary and conclusions. 

2. The Closest Point Method 

We now review the ideas behind the Closest Point Method [17] which are relevant to the 
problem of finding Laplace-Beltrami eigenvalues and eigenfunctions. 

The representation of the underlying surface is fundamental to any numerical method 
for PDEs on surfaces. The Closest Point Method relies on a closest point representation of 
the underlying surface. 

Definition 1 (Closest point function). Given a surface <S ; cp(sc) refers to a (possibly 
non-unique) point belonging to S which is closest to x. 

The closest point function, defined in a neighborhood of a surface, gives a representation 
of the surface. This closest point representation allows for general surfaces with boundaries 
and does not require the surface to have an inside/outside. The surface can be of any 
codimension [17], or even of mixed codimension [20]. 

The goal of the Closest Point Method is to replace a surface PDE by a related PDE 
in the embedding space which can be solved using finite difference, finite element or other 
standard methods. In the case of the Laplace-Beltrami eigenvalue problem, this approach 
relies on the following result, which states that the Laplace-Beltrami operator A^ may be 
replaced by the standard Laplacian A in the embedding space under certain conditions. 

Theorem 1. Let S be a smooth closed surface in R. d and u : S —> R be a smooth function. 
Assume the closest point function cp(x) is defined in a neighborhood <zM. d of S. Then 



Note that the right-hand side is well-defined because u(cp(-)) can be evaluated at points both 
on and off the surface. 

This result follows from the principles in [17]. 

Because the function u(cp(x)), known as the closest point extension of u ) is used through- 
out this paper, we make the following definition. 

Definition 2 (Closest point extension). Let S be a smooth surface in R d . The closest 
point extension of a function u : S — >> K. to a neighborhood Q of S is the function v : Q —> R 
defined by 



Asu(x) = A(u(cp(x))) for x G S. 



(2) 



v{x) = u(cp(x)). 



(3) 
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3. The embedded eigenfunction problem 

Our objective is to develop a simple, effective method for solving the following surface 
eigenvalue problem: 

Problem 1 (Laplace-Beltrami eigenvalue problem). Given a surface S , determine the 
eigenf unctions u : S M and eigenvalues A satisfying 

-A s (u(x)) = \u(x), for xeS. (4) 

If the surface S is open, then the problem will also have boundary conditions: we address 
this in Section 5. 

In this section, we assume that a closest point representation of the surface is available 
and consider two associated embedding problems. We will see that the first, a direct exten- 
sion of (2), leads to an ill-posed problem. However a straightforward modification of this 
problem leads to our second embedding problem, which we show is equivalent to (4). 

3.1. A first try 

We first consider the following embedded eigenvalue problem, which is directly motivated 
by (2). 

Problem 2 (Ill-posed embedded eigenvalue problem). Determine the eigenf unctions 
v : Q C M. d —> K. and eigenvalues A satisfying 

-A{v{cp{x))) = Xv{x), (5) 

in a neighborhood Q CM? of S. 

Solutions to Problem 1 and Problem 2 are closely related. Every solution to the embed- 
ding problem, restricted to the surface, is a solution to the surface problem. Conversely, 
except for the A = case, every surface eigenfunction corresponds to a unique solution of 
Problem 2. These results are established in Appendix A. 

Notably, the one-to-one correspondence between solutions breaks down for the A = 
case (the null-eigenspace) . Not only is this case significant in theory, in practice it makes 
Problem 2 ill-posed. 

3.1.1. The null-eigenspace 

The constant eigenfunction u{x) = c and A = is a solution to (4). Now consider a 
function on Q which agrees with u on the surface (but differs off the surface) 

v : Q —> R, such that v(x) = c for x e S, 

and note that v(x) is a null-eigenfunction of (5). Because v{x) is arbitrary for x G Q\S, the 
set of null-eigenfunctions for (5) is much larger than the set of null-eigenfunctions for (4). In 
fact, the set of null-eigenfunctions for (5) is infinite-dimensional: any (linearly independent) 
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change off the surface gives a new linearly independent eigenf unction. For points off of the 
surface, these null-eigenfunctions need not even be smooth. 

This example demonstrates the essential flaw of Problem 2: when A = 0, only the surface 
values of eigenfunctions are determined by (5) (i.e., eigenfunctions can take on arbitrary 
values elsewhere) . Not surprisingly, the infinite-dimensionality of the null-eigenspace causes 
problems for numerical methods based on approximating (5). For example, a large number 
of eigenfunctions with near-zero eigenvalues are observed in the numerical experiment of 
Figure 2 in Section 4.3 below. 

3.2. The fix: a modified embedded eigenvalue problem 

To avoid the null-eigenspace found in Problem 2, we consider a modified embedded 
eigenvalue problem. Our approach uses the split operator introduced in [20]: 

Definition 3 (Operator A*). Given fl C R d containing a surface S and a function v : 
Q —> R, the operator A* is defined as 

2d r 1 
Af v(x) := A(v(cp(x))) v(x) — v(cp(x)) . (6) 

where 0<£<1. 

The factor 2d (twice the dimension of the embedding space Q) is for later notational conve- 
nience. We can view Afv as A(v(cp)) plus a penalty for large change in the normal direction: 
specifically Af v will be large if \v(x) — v(cp(x))\ is not 0{e 2 ). Using this operator, we pose 
another embedded eigenvalue problem: 

Problem 3 (Regularized embedded eigenvalue problem). Determine all eigenfunc- 
tions : (] C K d 4 R and eigenvalues A satisfying 

-A*v(x) =\v(x), (7) 

in an embedding space ftcK d containing the surface S. 

For eigenvalues A < |f , we can show a one-to-one correspondence between Problem 1 
and Problem 3. 

Theorem 2 (Equivalence of two eigenvalue problems). Suppose S is a smooth sur- 
face embedded in and that Q C M d is a neighborhood of the surface. Then, for every 
eigenf unction u : C 2 {S) 4io/ (4) with eigenvalue A < ; there exists a unique eigenf unc- 
tion v : Q —> M of (7) with eigenvalue A which agrees with u on S. The eigenfunction v is 
given by 

A(u(cp(x))) + %u(cp(x)) 

< x ) = • (8) 



-A + 



2d 

^2 



Conversely, for every eigenfunction v{x) of (7) with eigenvalue \, the restriction of v{x) to 
S is an eigenfunction of (4) with eigenvalue A. 

PROOF. Existence and uniqueness of v follow directly from the condition that v agrees with 
u on S and the definition of Afv. The converse follows by choosing x G S in (7) and 
applying Theorem 1. 
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Remark 1. Note that when A = in (7), v is determined in Q by its values on the surface, 
avoiding the null-eigenspace problem of (5). 

Remark 2. If one tries to compute the eigenvalues A > |f of Problem 3, one encounters 
ill-posedness similar to the ill-posed embedded eigenvalue problem (5). This is confirmed 
by the numerical results in Section 6.1.3. The shift of the null-space by |f is due to the 
additional term |f (v(x) — v(cp(x))^j in Problem 3. As we shall see in the next section 
and Section 6.1.3, this ill-posedness will not be an issue for practical purposes since the 
parameter e can be chosen proportional to the grid spacing. 

Example 1. Let S be a circle of radius R embedded in 2D and consider the eigenf unction 
u = cos(^/\R9) with eigenvalue A. Applying Theorem 2, we get 



and we note that indeed v(x) = u(x) on the surface. 

3.2.1. Discretizing the regularized operator: choice of e 

In this work, we choose e = Ax, where Ax is the underlying grid spacing. In two 
dimensions, with this choice of s and a centered five-point discretization of the Laplacian 
operator, we begin discretizing the regularized operator (6) to obtain 

A*u(x, y) « — ^ [ - Au(x, y) + ?i(cp(x + Ax, y)) + u(cp((x - Ax, y)) + 



and analogously in higher dimensions. Note that by choosing e = Ax, the resulting di- 
agonal term (i.e., the "center" of the finite difference stencil) is without a closest point 
extension: where diagonal terms appear in the discretization, consistency does not require 
an interpolation to the closest point on the surface. See [20] for details. 

While (9) is an approximation to the regularized operator (6), it is not completely discrete 
since, for example, cp(x + Ax, y) is generally not a grid point. In Section 4.3 we will complete 
this discretization and solve the resulting eigenvalue problem. 

4. Numerical Discretization 

In this section, we discretize A* it following the approach given in [20]. The eigenvalues 
of the resulting matrix, computed using standard techniques, approximate the eigenvalues 



Following [20], our computational domain is a narrow band of m grid points L = 
{asi, #2, • • • , x m } enveloping the surface S. We approximate a function u at all points in L 
as the vector u G M m where U{ ~ u(xj) for each X{ G L. 




u(cp(x, y + Ax)) + ?i(cp(x, y — Ax))] , (9) 



of (4). 
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Figure 1: A discrete closest point exten- 
sion being applied to extend u(cp(x)) us- 
ing degree p = 3 (bicubic) interpolation. 
The shaded region indicates the interpo- 
lation stencil for cp(cc). 



4.1. Discrete closest point extensions 

To choose the list of points L, we consider the problem of evaluating u(cp(xj)) for a 
particular point Xj. As shown in Figure 1, cp(xj) is not typically a grid point, so we 
approximate the value of u(cp(xj)) using interpolation of the values of u at neighboring grid 
points. This interpolation has a certain stencil associated with it [20] and we choose L such 
that it contains the union of all the required interpolation stencils (an algorithm for doing 
this is given in [20]). 

In this work, we use Barycentric Lagrange polynomial interpolation [21] which is linear 
in the values of u. This allows us to express the closest point extension to any chosen set 
of points as a multiplication by an extension matrix E, where each row of E represents one 
discrete closest point extension. 

4-2. Discretizing the Laplace operator 

Combining the interpolation step, given by the matrix E, with a standard, symmetric 
finite difference discretization of A (e.g., the standard five-point stencil in 2D or the standard 
nine-point stencil in 3D), yields a discretization of A (u{cp(L))): 



A(u(cp(L))) « A h Eu = 



where M = A^E. is anmxm matrix intended to approximate the Laplace-Beltrami operator. 

However, the matrix M is based on discretizing A(^(cp(-))) and we know the latter 
cannot capture the eigenmodes of As correctly because of the issue of the null-eigenspace 
discussed in Section 3.1.1. Thus, it is not surprising to find that the discrete operator M also 
does a poor job of approximating the Laplace-Beltrami operator. For example, Figure 2 
shows that M has many eigenvalues close to zero, including some with imaginary parts and 
negative real parts. Notably, the matrix M is also ill-suited for time-dependent calculations; 
for example, [20, 19] show that implicit methods built on the matrix M have very strict 
stability time-stepping restrictions and are not competitive with explicit schemes. 

4-3. A modified discretization 

We will approximate the operator Af using a matrix M as described next. This approach 
will yield a convergent algorithm for the eigenfunctions and eigenvalues of the Laplace- 
Beltrami operator As- 
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Figure 2: Eigenvalues and histograms of their distribution for the matrices M (left) and 
M (right). The geometry is a unit semi-circle in 2D with homogeneous Dirichlet boundary 
conditions, E with bicubic interpolation (p = 3), and a mesh spacing of Ax = ^. Note the 
large number of eigenvalues near zero for M, whereas M correctly captures the spectrum of 
1,4,9,16,25,.... 



Approximating Af is straightforward and follows from the definition (6). When approx- 
imating As at the point Xi with a finite difference scheme, we map only the neighboring 
points Xj of the stencil back to their closest points cp(xj) and use U{ itself instead of u(cp(xi)) 
[20]. This special treatment of the diagonal elements (which corresponds to e — Ax) yields 
a new m x m matrix [20] 



M := diag(A /l ) + (A h - diag(A /l ))E = 



diagA^ Ah - diagA h 



Figure 2 gives the results of an experiment which contrasts the behavior of M and M. We 
find that the spectrum of M matches that of A^ on a semi-circular domain in M 2 , whereas 
the spectrum of M does not. 

Remark. The computational band L must be chosen so that it contains the interpolation 
stencil for cp(x) for all x needed in the calculation. Figure 3 shows the computational 
grid for an egg-shaped curve S illustrating the list of grid points L. An algorithm for the 
construction of the list L is given in [20] . 

4-4- Computing eigenvalues and eigenmodes 

Given the discrete operator M which approximates the Laplace-Beltrami operator on a 
surface 5, we compute a spectral decomposition 

M = QAQ-\ (10) 
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Figure 3: Computational grids for the Closest Point Method for an egg-shaped surface S 
with Ax = 0.125 (left) and Ax = 0.0625 (right). Second-order finite differences and degree 
p = 3 interpolation. 

where diag(A) are the eigenvalues of M and the columns of Q are the eigenvectors. These 
eigenvalues and eigenvectors are the respective approximations of the eigenvalues and eigen- 
functions of the Laplace Beltrami operator on S. 

J1..J1..I. Implementation 

The final step of our algorithm — computing the spectral decomposition of a matrix M — is 
a well-known problem in numerical linear algebra (e.g., [22, 23]). For example, in Matlab 
we can use the function eig() to compute the complete decomposition or the function eigs() to 
compute only part of the spectrum. Matlab's eig() calls LAPACK routines [24] and eigs() 
calls ARPACK [25] which makes use of the sparsity of the matrix M (in fact, it only requires 
a function which returns a matrix- vector product, although in this work we explicitly form 
the matrix M). Many of our calculations are performed in Python using SciPy [26] and 
NumPy [27] where we also make use of ARPACK via scipy.sparse.linalg.eigen.arpack. 

The eigenfunctions are returned as vectors over the list of points L. In practice, the inter- 
polation scheme from the closest point extension can be reused to evaluate the eigenfunctions 
at any desired points on the surface for plotting or other purposes. 

4.4.2. Degree of interpolation 

For the Closest Point Method to achieve the full order of accuracy of the underlying 
finite difference scheme (say q) , the degree of interpolation p must be chosen large enough so 
that the interpolant itself can be differentiated sufficiently accurately. For the second-order 
Laplace-Beltrami problems in this work, we need p > q + 1 [17, 20]. 
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Figure 4: The function cp. For the 
point x ) the closest point cp(x) lies 
in the interior of the surface and for 
such a point, cp(x) = cp(x). However, 
for x g) the closest point cp(a^) lies on s 
the boundary of the surface and cp(a^) ^ 
lies in the interior, with cp(a^) and x g 
roughly equidistant from cp(x g ). 

5. Boundary conditions 

When applied to an open surface, the Closest Point Method propagates boundary values 
into the embedding space along directions normal to the boundary, yielding homogeneous 
Neumann boundary conditions [17]. An analogous method for Dirichlet boundary conditions 
is similarly straightforward: instead of propagating out the interpolated values at boundary 
points the prescribed boundary conditions are propagated out [17]. While these methods 
are simple, they are only first-order accurate, which is lower-order than the discretizations of 
the Laplace^Beltrami operator discussed in this paper. Fortunately, a simple modification of 
the closest point function can be introduced to obtain a second-order accurate discretization 
for boundary conditions of Neumann or Dirichlet type. 

Assume a smooth surface, and consider the following modification of the closest point 
function, 

cp(sc) := cp(# + 2(cp(x) — as)) = cp(2cp(#) — x). (11) 

As is illustrated in Figure 4, whenever cp(x) is a point in the interior of the surface, the 
line between 2cp(x) — x and cp(x) is orthogonal to the surface. This implies that cp(x) = 
cp(2cp(#) — x) — cp(cc), at least in a neighborhood of the surface. Conversely, if cp(a^) ^ 
cp(a^) for a point x g then cp(a^) is an interior point of the surface 4 and cp(a^) is a boundary 
point. For a straight line or a planar surface, ^(cp(a^)) gives the mirror value for x g , while 
for a general, curved surface it gives an approximate mirror value. In correspondence to the 
terminology for codimension-zero regions with boundaries, we call a point x g a ghost point 
if cp(a^) 7^ cp(a^) (and note this terminology differs slightly from [20]). 

Thus, replacing cp(x) by cp(x) in the Closest Point Method does not change the treat- 
ment of interior points. At ghost points, however, (approximate) mirror values are used. 
This yields a second-order approximation of homogeneous Neumann conditions; no other 
modification of the method is required. Homogeneous Dirichlet boundary conditions are 
obtained by extending the function at ghost points by — u(cp(x g )). By analogy to second- 
order boundary conditions for codimension-zero regions, we observe that an approxima- 
tion of non-homogeneous Dirichlet conditions is obtained by extending at ghost points by 
u(xg) = 2u(cp(xg)) — u(cp(x g )), where u(cp(x g )) is di prescribed value on the boundary. 



4 At least for x g in a neighborhood of a sufficiently well-behaved surface: for example, for x g far from 
one boundary of a curve segment, cp(2cp(x^) — x g ) might be another boundary point instead of an interior 
point. 
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Figure 5: Numerical convergence studies for eight of the Laplace-Belt rami eigenvalues (X n = 
n 2 ) of a closed egg-shaped curve in 2D. The dashed reference lines have slope two (left) and 
four (right). Note second-order accuracy using second-order finite differences with degree 
p = 3 interpolation (left) and fourth-order accuracy using fourth-order finite differences with 
degree p = 5 interpolation (right). The lack of smoothness in the fourth-order results plot 
may be because the underlying spline curve itself is only C 2 smooth; results on a circle (not 
included) are smoother. 



We emphasize that the previous discussion relates to boundaries of the surface and the 
treatment of boundary conditions imposed thereat, not to the boundary of the narrow band 
of points (e.g., in Figure 3) where no artificial boundary condition need be applied [17, 20]. 

6. Numerical Results 

We present a series of numerical examples to demonstrate the effectiveness of our ap- 
proach. 

6.1. Numerical convergence studies 

We consider the egg-shaped curve in Figure 3 with arclength 2n as a test case. The 
Laplace-Beltrami eigenfunctions and eigenvalues in this case are the same as those of u" = 
\u, where u is 2tt periodic, namely u = cos (ns + f3) and A = n 2 for n e Z > where s 
represents arclength and f3 is a phase shift. The closest point function for this curve was 
determined using a numerical optimization procedure based on Newton's method. 

Figure 5 shows convergence studies in Ax for the first few eigenvalues on the egg-shaped 
domain. For second-order finite differences and degree-three interpolation we observe second- 
order convergence. That is, the Closest Point Method approximates the eigenvalues with 
error 0(Ax 2 ). Note that the error increases (and this is true even if one measures relative 
error) for the larger eigenvalues, but still shows second-order convergence. Figure 5 also 
shows that using fourth-order finite differences and degree-five interpolation, we get fourth- 
order accurate approximation to the eigenvalues. 
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(a) 2nd-order FD, p = 3 interp. (b) 4th-order FD, p = 5 interp. 

Figure 6: A curve and the corresponding computational grids for the Closest Point Method 
with second-order (left) and fourth-order (right) finite difference schemes with Ax = 0.125. 
Grid points marked with small crosses are the ghost points involved in implementing bound- 
ary conditions. 



In further numerical tests (not included), we found that degree p = 2 interpolation 
with second-order finite differences, and degree p = 4 interpolation with fourth-order finite 
differences, give approximately second- and fourth-order convergence, respectively. These 
are better than expected by one order of accuracy, as originally observed in [20] . 

6.1.1. Boundary Conditions 

In this section we verify using convergence studies that our treatment of boundary con- 
ditions introduced in Section 5 achieves the expected order of accuracy. For our tests, we 
use the curve shown in Figure 6 parameterized as (x(t),y(t)) = (t, cost) for t G [1/4,4]. 
We apply homogeneous Neumann and Dirichlet boundary conditions to the ends. Again, 
the exact eigenvalues and eigenfunctions can be determined analytically by considering the 
corresponding problem on an interval. For the former, the arclength of the curve is required 
and can be determined in terms of elliptic integrals. 

The original Closest Point Method, which does not explicitly treat boundaries, gives 
a first-order approximation to homogeneous boundary conditions [17]. See, for example, 
Figure 7 where it is observed that this trivial treatment of the boundaries limits the overall 
accuracy in the eigenvalues to first-order. Figure 7 also gives results using the new a cp" 
approach for Neumann boundary conditions described in Section 5; we find that this method 
attains the second-order accuracy of the underlying Cartesian finite difference scheme. Note 
that the cp approach is also effective for Dirichlet boundary conditions. For example, in 
Figure 8 we find second-order results for homogeneous Dirichlet boundary conditions using 
—u(cp(x g ) for ghost points (as in Section 5). 

The cp approach is designed to give second-order approximations to boundary conditions. 
Thus it is expected, and observed in Figure 8, that even with higher-order finite difference 
schemes and high-degree interpolation, the results will generally be second-order accurate in 
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Figure 7: Numerical convergence studies for the first eight non-zero Laplace-Beltrami eigen- 
values of a curve in 2D with Neumann boundary conditions applied at both ends. The left 
figure uses the original closest point function to impose a homogeneous Neumann boundary 
condition, and we note the results are only first-order accuracy (dashed reference line has 
slope one). The right figure uses the modified "cp" approach and exhibits second-order 
accuracy (dashed reference line has slope two). Both computations use second-order finite 
differences and degree p = 3 interpolation. 



Ax 



size of M 



tf(M) 



0.25 
0.125 
0.0625 
0.03125 
0.015625 
0.0078125 



76 x 76 
140 x 140 
268 x 268 
524 x 524 
1036 x 1036 
2060 x 2060 



289 

1154 

4608 

19304 

75543 

326633 



Table 1: Condition numbers (2- 
norm) for the matrix M. Tested 
for the Laplace-Beltrami operator 
on a open curve embedded in 2D 
with Dirichlet boundary conditions at 
both ends, using second-order finite 
differences, degree p = 3 interpola- 
tion, and the "cp" treatment of the 
boundary conditions. 



the presence of boundary conditions. Third- and higher-order approximation of boundary 
conditions may also be contemplated; while we do not investigate such methods here we 
note that approximations of this type will require a replacement for cp that incorporates 
the curvature of the surface near the boundary. 

6.1.2. Conditioning 

Table 1 shows that the condition number of the matrices used in our computations scales 
like 0(^3). Thus the conditioning is the same as for standard Cartesian finite difference 
schemes. The "cp" treatment of boundary conditions does not have a significant effect on 
conditioning. 

6.1.3. Behavior for high-frequency modes 

The large number of spurious eigenvalues near i n Figure 9 reflect the singular be- 
havior of Equation (8) at A = (c.f., Figure 2 where the problem instead occurs near the 
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Figure 8: Numerical convergence study for the first few eigenvalues on a curve with Dirichlet 
boundary conditions applied at both ends. The boundaries are dealt with using the a cp" 
and — u discretization with degree p interpolation. The left figure uses second-order finite 
differences (and degree p = 3 interpolation) and achieves second-order accuracy (the dashed 
reference lines have slope two). The right figure uses fourth-order finite differences (and 
degree p = 5 interpolation) but the overall accuracy is limited by the second-order treatment 
of the boundary conditions. 



origin). Notably, it is these spurious eigenvalues which possess nonzero imaginary compo- 
nents, reflecting the fact that the eigenvalue problem (7) is not self-adjoint. Because these 
spurious eigenvalues correspond to highly oscillatory eigenfunctions which are quite close to 
the Nyquist frequency (i.e., with eigenvalues near ^2), it is not surprising that numerical 
difficulties arise for such modes. Nonetheless, any particular higher-frequency modes may 
be resolved by appropriately refining Ax. 

6.2. Eigenf unction computations 

The following computations were performed in Matlab and SciPy [26]. Visualizations 
were carried out with Matlab, VTK [28] and MayaVi [29]. 

6.2.1. Hemispherical harmonics 

Figure 10 shows several Laplace-Beltrami eigenmodes of a unit hemisphere with homo- 
geneous Neumann boundary conditions on the equator. 

6.2.2. Eigenkaninchen 

Figure 11 shows Laplace-Beltrami eigenmodes of the surface of the Stanford Bunny [30]. 
The closest point function for this and other triangulated surfaces can be computed in an 
efficient, straightforward manner [18]. 

6.2.3. Mobius strip 

Figure 12 shows some Laplace-Beltrami eigenmodes of a Mobius strip. Dirichlet bound- 
ary conditions are applied at the boundary of the Mobius strip. The surface is embedded in 
M 3 and the closest point function for each grid point is computed from a parameterization 
using a numerical optimization procedure minimizing the square of the distance function. 
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256 



384 



-3 







512 



1024 



1536 



-15 LL 




2048 



4096 



6144 



real(A) 



real(A) 



real(A) 



Figure 9: Computed spectra of the Laplace-Beltrami operator on a unit circle using Ax = 
|, and ^ (from left-to-right). The method uses second-order finite differences and degree 
p = 3 interpolation. Note a large number of extraneous complex eigenvalues near (i-e., 
256, 1024, and 4096): these correspond to the singularity in (8) and can be controlled by 
further resolving Ax. 



Figure 10: The left image shows several hemispherical harmonics with eigenvalue 20 com- 
puted with second-order finite differences, degree p = 3 interpolation, and Ax = ^. The 
right image is a convergence study for the first six eigenvalues (A n = n(n + 1)), also using 
second-order finite differences and degree p = 3 interpolation (the dashed reference line has 
slope two). 




l/Ax 
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Figure 11: A selection of eigenkaninchen: Laplace-Beltrami eigenmodes on the surface of 
the Stanford Bunny [30]. The results are computed via the Closest Point Method using 
second-order finite differences, degree p = 3 interpolation, and Ax = 0.1, where the bunny 
is roughly two units long. 



Figure 12: A selection of Laplace-Beltrami 
eigenmodes of a Mobius strip computed with 
the Closest Point Method. Calculations use 
second-order finite differences, p — 3, Ax = 
0.1, and the Mobius strip is about 2 units in 
"diameter" . 




Note that the non-orientable nature of the Mobius strip poses no difficultly for the Closest 
Point Method. 



6.2.4- L- shaped Domain 

The Closest Point Method works on surfaces of various codimension [17, 20] and indeed 
solid shapes in K 2 or M 3 are surfaces of codimension-zero. Figure 13 shows an eigenmode 
computation on an L-shaped domain, where zero Dirichlet boundary conditions are imposed 
using the "cp" technique described in Section 5. In the interior of a solid, cp(x) = x and so 
no interpolations are needed. Furthermore, for a grid-aligned L-shaped domain, for any X{ 
outside the domain, cp(a^) turns out to be another grid point Xj (located inside the domain 
as if the perimeter were a mirror) so no interpolation step is necessary. Thus in this special 
"corner case" , the Closest Point Method reduces to a standard finite difference computation 
using ghost points to mirror the values along the perimeter. Interestingly, these reductions 
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WTO*. 



Figure 13: The first 15 Laplacian eigen- 
modes of an L-shaped domain with 
Dirichlet boundary conditions. The re- 
sults are computed using the Closest 
Point Method with second-order finite 
differences using Ax = 1/40 (the domain 
is two units wide). 



*T tTvitI f 




Figure 14: First ten eigenmodes of continental Africa. The Earth has unit diameter and 
the computation uses second-order finite differences with Ax = 1/40 and degree p = 3 
interpolation. 



happen automatically: no change in the code is needed. 

6.2.5. Continental Africa 

Some Laplace-Beltrami eigenmodes of continental Africa are shown in Figure 14. The 
results are computed directly on the surface of the Earth (assumed to be a sphere). Homo- 
geneous Dirichlet boundary conditions are applied to the coastline. Finding the closest point 
function involves projecting onto a bitmapped image of the Earth [31] where the continent 
was first manually segmented. It is interesting to note that these eigenmodes match very 
closely the first ten eigenmodes of the L-shaped domain in Figure 13. 



7. Summary and Conclusions 

Through a series of convergence studies and computational examples, we have shown 
that the Closest Point Method is an effective method for computing the spectra and eigen- 
functions of the Laplace-Beltrami operator on rather general surfaces. The basis of our 
approach is the embedded eigenvalue problem (7) which, when discretized using standard, 
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centered finite difference methods and Lagrange interpolation, yields a nonsymmetric ma- 
trix eigenvalue problem which can be solved using standard software. Fortunately, this lack 
of symmetry is not a concern in many practical situations since only the highest frequency 
modes have a significant imaginary component. We are currently investigating a finite ele- 
ment Closest Point Method which would lead to symmetric matrices. 

For eigenvalue problems on open surfaces with Dirichlet or homogeneous Neumann 
boundary conditions, we have introduced second-order accurate schemes. These methods 
require only a simple change to the closest point extension and are straightforward to im- 
plement. 

Although we have focused on the Laplace-Beltrami operator in this work, the Closest 
Point Method is applicable to many other surface differential operators [17, 18, 19, 20]. This 
suggests that the approach presented here may be applicable to a larger class of eigenvalue 
problems. 

Finally, the techniques developed here are quite general and are applicable beyond surface 
eigenvalue problems. For example, future work will include solving surface Poisson problems 
of the form —Agu f. 
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Appendix A. Additional theorems 

As mentioned in Section 3.1, every solution to the ill-posed embedding problem (5), 
restricted to the surface, is a solution to the surface eigenvalue problem. Conversely, except 
for the A = case, every surface eigenfunction corresponds to a unique solution of the 
embedding problem. These results are established by the following theorems. 

Theorem 3. Suppose v{x) and A are a solution to the embedding eigenvalue problem (5) 
and S is a smooth surface. Then u : S —> M defined by u(x) = v(x) for x e S is an 
eigenfunction of (4) with eigenvalue A. 

PROOF. This is a direct consequence of Theorem 1. 

Theorem 4. Let S be a smooth surface. For every non-constant solution u and A of (4), 
there exists a unique (up to a multiplicative constant) eigenfunction v of (5) with eigenvalue 
A which agrees with u on S. This eigenfunction is given by v(x) = — jA(u(cp(x))). 

Proof. Using our hypothesis that v agrees with u on <S, we may solve for v in (5) to obtain 
the result. 
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